/*******************************************************************************

This code file produces Figure 2, "Zeroth and First Stages of Simulated Instrument GMM."

*******************************************************************************/

*** Manage settings

	run "$dir/code/modules/settings.do"

********************************************************************************
* Simulated instruments GMM
********************************************************************************
	
*** Collapse to means of matched-building features and merge

	use "$data/raw/matched_buildings.dta", clear
	
	format bbl %18.0f
	
	gen mvalue_pu = assesstot/(unitsres*frac_assess*underassess)
	
	gen sh_exempt = exempttot/assesstot
	
	gcollapse (mean) unitsres bldgarea numfloors dtaxrate_onsite builtfar condo mvalue_pu sh_exempt, by(bbl)
	
	foreach v of varlist unitsres bldgarea numfloors dtaxrate_onsite builtfar condo mvalue_pu sh_exempt {
		rename `v' `v'_sim
	}
	
	merge 1:1 bbl using "$data/clean/cleaned_data.dta", nogen
	
	gen mvalue_pu = assesstot/(unitsres*frac_assess*underassess)
	
*** Regressions

	foreach v of varlist unitsres bldgarea builtfar mvalue_pu condo {
	
		quietly summ `v', d
		local p1 = r(p1)
		local p99 = r(p99)
	
		reg `v' `v'_sim if `v'_sim > `p1' & `v'_sim < `p99', r cl(nta)
		store_est_tpl using analysis_t11.csv, coef(`v'_sim) name(`v') all
		
	}
	
*** Binscatters

	** Zeroth stage: Are the buildings similar on observables?
	
	replace bldgarea = bldgarea / 1000
	replace bldgarea_sim = bldgarea_sim / 1000
	
	replace mvalue_pu = mvalue_pu / 1000000
	replace mvalue_pu_sim = mvalue_pu_sim / 1000000
	
	foreach v of varlist unitsres bldgarea builtfar mvalue_pu condo {
	
		quietly summ `v'_sim, d
		local p1 = r(p1)
		local p99 = r(p99)
		
		capture gr drop gr_`v'
		
		if "`v'" == "unitsres" {
			
			binscatter `v'_sim `v' `v'_sim if `v'_sim > `p1' & `v'_sim < `p99', ///
				n(100) line(lfit) lcolor(gs9 none) mcolor(none black) ///
				aspect(1) legend(off) name(gr_`v') ylabel(, nogrid) graphregion(lcolor(white)) ///
				ytitle("Num. of Residential Units, Actual") xtitle("Num. of Residential Units, Matched")
			
		}
		
		if "`v'" == "bldgarea" {
			
			binscatter `v'_sim `v' `v'_sim  if `v'_sim > `p1' & `v'_sim < `p99', ///
				n(100) line(lfit) lcolor(gs9 none) mcolor(none black) ///
				aspect(1) legend(off) name(gr_`v') ylabel(, nogrid) graphregion(lcolor(white)) ///
				ytitle("Building Area (K Sq. Ft.), Actual") xtitle("Building Area (K Sq. Ft.), Matched")
			
		}
		
		if "`v'" == "condo" {
			
			binscatter `v'_sim `v' `v'_sim if `v'_sim > `p1' & `v'_sim < `p99', ///
				n(100) line(lfit) lcolor(gs9 none) mcolor(none black) yscale(range(0 1)) ylabel(0(0.2)1) ///
				aspect(1) legend(off) name(gr_`v') ylabel(, nogrid) graphregion(lcolor(white)) ///
				ytitle("Condominium Share, Actual") xtitle("Condominium Share, Matched")
			
		}
		
		if "`v'" == "builtfar" {
			
			binscatter `v'_sim `v' `v'_sim  if `v'_sim > `p1' & `v'_sim < `p99', ///
				n(100) line(lfit) lcolor(gs9 none) mcolor(none black) ///
				aspect(1) legend(off) name(gr_`v') ylabel(, nogrid) graphregion(lcolor(white)) ///
				ytitle("Floor Area Ratio, Actual") xtitle("Floor Area Ratio, Matched")
			
		}
		
		if "`v'" == "mvalue_pu" {
			
			binscatter `v'_sim `v' `v'_sim if `v'_sim > `p1' & `v'_sim < `p99', ///
				n(100) line(lfit) lcolor(gs9 none) mcolor(none black) xscale(range(0 1)) xlabel(0(1)4) ///
				aspect(1) legend(off) name(gr_`v') ylabel(, nogrid) graphregion(lcolor(white)) ///
				ytitle("Value Per Unit ({c $|}M), Actual") xtitle("Value Per Unit ({c $|}M), Matched")
			
		}
		
	}
	
	summ dtaxrate_onsite_sim, d
	local p1 = r(p1)
	local p99 = r(p99)
	
	preserve
	egen grp = cut(dtaxrate_onsite_sim), group(100)
	
	collapse dtaxrate_onsite dtaxrate_onsite_sim if dtaxrate_onsite_sim > `p1' & dtaxrate_onsite_sim < `p99', by(grp)
		
	capture gr drop gr_dtaxrate
	tw scatter dtaxrate_onsite dtaxrate_onsite_sim, mcolor(black) ///
		|| lfit dtaxrate_onsite_sim dtaxrate_onsite_sim, lcolor(gs9) ///
		ylabel(, nogrid) graphregion(color(white)) range(0 0.8) xscale(range(0 0.8)) xlabel(0(0.2)0.8) ///
		aspect(1) legend(off) name(gr_dtaxrate) ytitle("421-a Incentive, Actual") xtitle("421-a Incentive, Matched")
	restore

	gr combine gr_unitsres gr_bldgarea gr_builtfar gr_condo gr_mvalue_pu gr_dtaxrate, ///
		graphregion(color(white))
	
	graph export "$figs/siminst_zeroth_stage.pdf", replace
	graph export "$figs_overleaf/siminst_zeroth_stage.pdf", replace
